Evidence of spin liquid with hard-core bosons in a square lattice 



(N 

o 

(N 

< 
o 



C/3 

: 

wo. 

^— > ■ 
c : 

cr 

+-» ! 



i 

c 

o 
o 



> 

00 

in 

\o 

O 



X 



Y.-H. Chan 1 ' 2 and L.-M. Duan 1 ' 2 
1 Department of Physics, University of Michigan, Ann Arbor, Michigan 48109, USA and 
2 Center for Quantum Information, HIS, Tsinghua University, Beijing, China 

We show that laser assisted hopping of hard core bosons in a square optical lattice can be described 
by an antiferromagnetic J1-J2 XY model with tunable ratio of J2/J1 ■ We numerically investigate the 
phase diagram of the J\-J% XY model using both the tensor network algorithm for infinite systems 
and the exact diagonalization for small clusters and find strong evidence that in the intermediate 
region around J2/J1 ~ 0.5, there is a spin liquid phase with vanishing magnetization and valence 
bond orders, which interconnects the Neel state on the J2 <C Ji side and the stripe antiferromagnetic 
phase on the J2 3> Ji side. This finding opens up the possibility of studying the exotic spin liquid 
phase in a realistic experimental system using ultracold atoms in an optical lattice. 

PACS numbers: 



A spin liquid phase is an exotic state of matter that 
does not break any symmetry of the Hamiltonian and 
has no conventional order even at zero temperature [ij. 
A number of microscopic Hamiltonians with frustrated 
quantum magnetic interaction could support a spin liquid 
phase 1-6]. In particular, very recently, numerical inves- 
tigations based on complementary methods have found 
strong evidence that the antiferromagnetic J1-J2 Heisen- 
berg model may have a spin liquid phase in a square 
lattice 0, Q< On the experimental side, several mate- 
rials are suspected to be in a spin liquid phase at very 
low temperature [l|. However, due to complication of 
physics in these materials, it is hard to make a direct 
connection of the prediction from the simplified micro- 
scopic models and the phenomenology observed in real 
materials Ultracold atoms in an optical lattice pro- 
vides a clean platform to realize microscopic models to 
allow for controlled comparison between theory and ex- 
periments 0, HJ. Proposals have been made to imple- 
ment the frustrated magnetic models in an optical lattice 
0, [13] and various required configurations of the optical 
lattices have been realized experimentally ■ However, 
the direct magnetic Heisenberg coupling, which comes 
from the higher-order super-exchange interaction, is very 
weak under typical experimental conditions 0, [ll| . It is 
still very challenging to reach the extremely low tempera- 
ture required to observe the ground state of the magnetic 
Heisenberg model in an optical lattice. 

In this paper, we show strong evidence that a spin liq- 
uid phase can emerge in an antiferromagnetic J\-Ji XY 
model in a square lattice. The calculations are based on 
two complementary methods: the recently developed ten- 
sor network algorithm applied directly to infinite systems 
12 , [HI and the exact diagonalization of small clusters 
which is combined with the finite size scaling to infer the 
phase diagram [l4| ■ Both methods suggest that in a small 
region around J2/J1 ~ 0.5, magnetization and valence 
bond solid orders all vanish, indicting a spin liquid phase 
as the ground state. Different from a Heisenberg model, 
a XY model can be realized with hard-core bosons in 
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FIG. 1: (Color online) (a) Implementation of the Ji — Ji XY 
model with cold bosons in a bi-partite square optical lattice, 
where the J2 coupling is due to the atomic hopping in the 
same sub-lattice, and the Ji coupling is induced by the three 
Raman laser beams (the direct Ji hopping of the atoms is 
turned off by the large potential shift between the two sub- 
lattices), (b) The configuration of the wave- vectors for the 
three Raman laser beams. 



an optical lattice. Through control of the laser assisted 
hopping in a square lattice (lsj . we propose a scheme 
to implement the effective antiferromagnetic couplings 
for both the neighboring and the next neighboring sites 
with a tunable ratio of J%j J\. In this implementation, 
both J2 and J\ are determined by the hopping rates of 
the hard-core bosons in an optical lattice, which is much 
larger than the conventional super-exchange interaction 
for ultracold atoms in the Heisenberg model 0, [llj . The 
large J1-J2 couplings open up the possibility to exper- 
imentally realize this model and observe its spin liquid 
phase based on the state-of-the-art technology. 

The J1-J2 XY model is represented by the Hamilto- 



H 



Ji y^XjXj 

(id) 



J2 (XiXj 

((i,D) 



■YiYi), (1) 



where X, Y represent the Pauli operators a x and a y 



J) 



and {{i, j)) denote respectively the neighboring and the 
next neighboring sites in a square lattice as shown in Fig. 
1(a). To realize this model with hard core bosons, we 
consider ultracold atoms in different hypcrfinc spins \a) 
and j b) loaded into alternating square lattices A and B as 



2 



shown in Fig. 1. This configuration can be experimen- 
tally realized with the spin-dependent lattice potential 
(l6| . Atoms in spins \a) (or \b)) freely tunnel in the lat- 
tice A (or B) with the hopping rate t, however, a direct 
hopping between the A, B lattices is forbidden due to the 
spin-dependent potential shift. Instead, the inter-lattice 
hopping is introduced by the laser induced Raman tran- 
sition as shown in Fig. 1(a). We use three Raman beams, 
with wave- vectors ki, k2, and k3 and Rabi frequencies 
fil, Oa, and r?3, respectively. The directions of the laser 
beams are shown in Fig. 1(b) with Aki2 = ki— k2 = k^y 
and Aki3 = ki — k3 = kpjb. The lase induced inter- 
lattice hopping rates for the neighboring sites are then 
given by t x = J w* (xi,yi) Sltto 3 e ikAa >w {x l+ i,y l ) dxdy, 
and t 9 = J w* (xi,yi)rt\Q,2e lkAV w (x u y t+1 ) dxdy, for 
the hopping along the x, y directions, respectively. As- 
sume H3 = — O2 and the Wannier function w{xi,yi) 
symmetric along the x, y directions, we have t x = —t y = 
t' (we can always choose t' > by setting an appropriate 
relative phase between Jli and H3). If the on-site atomic 
repulsion U satisfies U 3> t, t', we have the hard-core con- 
straint with at most one boson per site. The hard-core 
bosons in this square lattice are then described by the 
Hamiltonian 



H = t'J2 4b, -f ]T a\bj-t 



E 



(4aj + blb j ) + H.c. 

(2) 

The hard core bosons <2j , bj satisfy the same commuta- 
tors as the Pauli operators a~ , a~ , so with the mapping 



a-i — i 
and a 



tj- and b-j 



for the odd numbers of rows, 
for the even numbers 



and bj — > — erj 
of rows, the Hamiltonian (2) is mapped to the J1-J2 XY 
model in Eq. (1) with J x = t'/2 > and J 2 = t/2 > 0. 
Apparently, the ratio J2/J1 is tunable by changing the 
magnitude of the Rabi frequencies H|$l3. 

In the following, we calculate the phase diagram of the 
Hamiltonian (1) as a function of the dimensionless pa- 
rameter J2/J1 (Ji is taken as the energy unit). In the 
limit J2/J1 "C 1, the J\ term dominates and the ground 
state is magnetized with a Neel order at the momen- 
tum k — (jr, 7r). In the opposite limit J2/J1 S> 1, the 
ground state has a stripe magnetic order at the momen- 
tum (71", 0) or (0, 7r), which minimizes the energy of the J2 
term. In the intermediate region with J2/J1 ~ 0.5, the 
Hamiltonian is highly frustrated with competing interac- 
tion terms. Our main purpose is to find out the phase 
diagram in this region through controlled numerical sim- 
ulations. 

Our numerical simulations are based on two compli- 
mentary methods: exact diagonalization (ED) for small 
clusters 
systems 



14 1 and tensor network simulation for infinite 
HE1. The ED method is limited by the clus- 



ter size, and we use extrapolation based on the finite-size 
scaling to infer the phase diagram for the infinite system. 
The tensor network algorithm is an recently developed 



simulation method inspired by quantum information the- 
ory [l2j]- It can be considered as an extension of the 
density matrix renormalization group (DMRG) method 
to the two dimensional case, replacing the matrix prod- 
uct state in the DMRG method with the tensor network 
state that better matches the geometry of the underly- 
ing lattice [12J ■ We use a particular version of the tensor 
network algorithms, the infinite projected entangled pair 
states (iPEPS) method [HI], which applies directly to in- 
finite systems using the translational symmetry To take 
into account the ordered states for the Hamiltonian (1) 
that spontaneously break the translational symmetry, in 
our simulation we take a unit cell (typically 2x2 and 
4x4) that is large enough to incorporate the relevant 
symmetry breaking orders 13] . We apply imaginary time 
evolution to reach the ground state of the Hamiltonian. 
To avoid being stuck in a metastable state, we take a 
number of random initial states for the imaginary time 
evolution and pick up the ground state as the one which 
has the minimum energy over all the trials. The accuracy 
of the iPEPS simulation depends on the internal dimen- 
sion D of the tensor network state. The simulation time 
scales up very rapidly with the dimension D, which lim- 
its D to a small value in practice. We typically take D 
between 4 to 6 in our simulation. 

Figure 2 shows the major result from the iPEPS 
simulation. First, we look at the average magnetiza- 
tion m s = (1/N S ) VXi 2 + Yi 2 + Zi 2 as a function of 
J2/J1, where the average is taken over the N s sites in the 
unit cell. The calculation shows that for small or large 
J2/J1, the ground states are magnetic (with the Neel or 
the stripe order, respectively), which is consistent with 
our intuitive picture. In the intermediate region with 
0.46 ^ J2/J1 ^ 0.54, there is a sudden drop of all the 
magnetic orders to a tiny value. Although the iPEPS 
method under a small dimension D could be biased to- 
ward a less entangled state, which is typically an ordered 
state, it would not be baised toward a disordered spin 
liquid state. So, when we see a big sudden drop of the 
magnetic orders from the simulation, it must be a real 
effect, strongly indicating there is a new phase in the in- 
termediate region with vanishing magnetic orders. The 
remaining small m s may be due to the finite dimension 
D and should vanish when D is scaled up. 

To figure out the property of the phase in the interme- 
diate region, we further check different kinds of valence 
bond solid orders. We calculate all the neighboring va- 
lence bonds (<7i • <jj) in the unit cell and the result is 
shown in Fig. 2. For a valence bond solid state, the 
spatial symmetry should be spontaneously broken for the 
valence-bond distribution. Figure 2 shows that in the en- 
tire region of J2/ J\, the valence bond distribution has the 
same symmetry as the underlying Hamiltonian, which in- 
dicates that the ground state of the Hamiltonian (1) has 
no valence bond solid orders. Together with the above 
calculation of the magnetic orders, this suggests that the 
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FIG. 2: (Color online) Average magnetization m s as a func- 
tion of 3-2.1 Ji . The insets show the spin configuration and the 
valence bond distribution (cr* ■ Oj) at J2/J1 = 0, 0.5, and 0.9 
obtained with the iPEPS on a 4 x 4 unit cell with D = 6. 
The width and color of the bonds are scaled such that the 
negative energy is shown by thicker bond with darker color 
and the positive energy is shown by thinner bond with lighter 
color and the length of the spin is proportional to its magnetic 
moment m a . 
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FIG. 3: (Color online) (a) Spin-spin correlation (<7i ■ cr,-) as a 
function of distance d along the diagonal direction at J%j J\ = 
0.1 (cross), 0.5 ( circle) and 0.9 (open diamond), (b) Semi-log 
plot of spin-spin correlation (Act* • A<7j) after subtracting the 
local averages, (c) Semi-log plot of dimer-dimer correlation 
(ADf AD") (a = x,y) as a function of distance d along the 
diagonal direction at J2/J1 = 0.5. 



Hamiltonian (1) has a spin liquid phase with no orders in 
the intermediate region with 0.46 ^ J2/J1 ^ 0.54. This 
spin liquid phase seems to have the same feature as the 
Z2 spin liquid in the intermediate coupling region of the 
J1-J2 Heisenberg model found in the recent numerical 
simulation 0, @. 

To further confirm this picture, we calculate the long- 
range spin correlation and dimer correlation with the 
iPEPS method and the result is shown in Fig. 3 for 
J2/J1 = 0.1,0.5 and 0.9. The spin correlation (<7j • Cj) 
is calculated along the diagonal direction x + y . Both 
the Neel and the stripe phases have long-range correla- 
tions, with constant or staggered values along the diag- 
onal direction. The intermediate phase has an exponen- 
tially decaying spin-spin correlation, which is in agree- 
ment with the behavior of the Z2 spin liquid phase with 
a finite spin gap [H, Hjj- The dimer operator Df is de- 
fined by Df — at ■ <Ji-\- a for the bond + a), where 
a = x or y denote the orientation of the dimer. In Fig. 
3(c), we show the dimer-dimer correlations (ADfADj) 
and (A£>f A£>p at J2/J1 — 0.5 along the diagonal direc- 
tion. The correlations are exponentially decaying with 
distance, in agreement with a spin liquid phase with no 
dimer orders. 

In the following, we present study of the Hamiltonian 
(1) with the complementary ED method, which provides 
further evidence for a spin liquid phase in the interme- 
diate region. To be consistent with the periodic bound- 
ary condition required for the finite size scaling and to 



incorporate the momentum k — {tt, tt) responsible for 
the Neel order, the size of the clusters for the ED is 
taken to 16, 20 and 32 sites. From the spin correla- 
tion (<Ti ■ (Tj), we calculate the corresponding static struc- 
ture factor m 2 s (k, N) = {1/N) e ik '( r '- r ^ (Act, • Act,-), 
where N is the size of the cluster and Ac; = a* — {(Ji). 
The Neel order and the stripe order correspond to peaks 
at k = {tt, tt) and {tt, 0), respectively. Finite-size clusters 
always have non-zero order parameters, and one needs 
to do finite size scaling, with a simple scaling formula 
m 2 s {k,N) — mf(k, 00) + a/y/N {VN corresponds to the 
linear size), to infer the value of m 2 {\s., 00) for the infi- 
nite system. In Fig. 4, we show the finite size scaling 
for m 2 s {k,N) at J2/J1 = 0,0.5, and 0.9 in three differ- 
ent regions. The results are consistent with the find- 
ings from iPEPS method, i.e., there is a stripe order 
with k = {tt, 0) at J2/J1 — 0.9 and a Neel order with 
k = (7r,7r) at J 2 /Ji = 0. At J 2 /Ji = 0.5, the finite- 
size scaling indicates a vanishing stripe order. However, 
at k = (7r,7r), the data become non-monotonic with N 
due to the shape of the cluster and the finite-size scaling 
becomes inconclusive in this case. The non-monotonic 
shape effect has also been observed in ED for the J1-J2 
Heisenberg model (l4| . 

To check for possible valence bond solid orders 
from ED, we similarly calculate the structure fac- 
tors m 2 {k,N) = {1/N) Y^ij e ik '( ri ~ r ^(ADfA.DJ) and 
m 2 p {k,N) = (l/A^)^, y .e ik -( r --^)((Ai 3 i AP J ), corre- 
sponding respectively to the dimer order Df and the pla- 
quette order Pj = {Qi + QJ )/2, where Qi {Q^ 1 ) is the 
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FIG. 4: (Color online) Finite size scaling of the magnetic order 
parameter at (a) k = (7r,7r) and (b) k = (ir, 0) at J2/J1 = 
(dot), 0.5 (square), and 0.9 (diamond). 
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FIG. 5: (Color online) Finite size scaling of (a) the dimer 
order parameter at k = (71-, 0) and (b) the plaquette order 
parameter at k = (7r,7r) at J2/J1 = (dot), 0.5 (square), and 
0.9 (diamond). 

clockwise (anticlockwise) cyclic permutation operator on 
the plaquette i with its explicit (lengthy) expression given 
in [171, [l8| . The rotational symmetry is always preserved 
at finite size, so we only need to check one component 
of the dimer order, say Df. At finite size, the structure 
factors peak at k = (71", 0) for the dimer order Df and at 
k = (tt, 7r) for the plaquette order Pj, however, an extrap- 
olation to the infinite system at these momenta as shown 
in Fig. 5 indicates vanishing dimer and plaquette orders 
in all three regions of J2/J1. This result, again, is in 
agreement with the finding from the iPEPS calculation. 

Before concluding the paper, we briefly discuss the ex- 
perimental signature of the three different phases for the 
Hamiltonian (1) in the implementation with hard-core 
bosons. The Neel ordered state and the strip phase cor- 



respond to Bose-Einstein condensates at the momenta 
k = (tt, 7r) and k = (tt, 0), respectively. The standard 
time-of-flight imaging measurement can then reveal the 
condensate peak at these nontrivial momentum points 
!8!j. The spin liquid phase, on the other hand, would not 
show any condensation peaks due to lack of magnetic 
orders. Furthermore, it has a spin gap which implies a 
charge gap in implementation with hard-core bosons. We 
therefore expect to see an incompressible phase at half 
filling, which is different from the Mott insulator state 
at the integer filling. It is also be distinguished from a 
charge density wave state since the density distribution 
in this case is still homogeneous without any solid order. 

In summary, we have proposed an experimentally fea- 
sible scheme to implement the J1-J2 XY model with 
ultracold hard-core bosons in a square optical lattice. 
Through detailed numerical simulation of this model us- 
ing two complementary methods, we find strong evidence 
that this model has a spin liquid phase in the interme- 
diate region of Jz/Ji- The proposed experimental im- 
plementation, with tunable ratio of Jij J\, opens up a 
realistic possibility to look for the long-pursued spin liq- 
uid phase in a well controlled Hamiltonian model. 
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